// STL
#include <tuple>                    // std::tie()
#include <utility>                  // std::pair<>, ::make_pair()
#include <vector>                   // std::vector<>

// stats++
#include "statsxx/distribution.hpp" // distribution::CDF


//
// DESC: Estimate the probability density function (PDF) of the probability distribution from which a set of data was drawn.
//
// -----
//
// NOTE: *CAUTION* Using finite-differences for a CDF from discrete data does *NOT* work well.
//
inline std::pair<
                 std::vector<double>, // PDF_x
                 std::vector<double>  // PDF_y
                 > distribution::PDF(
                                     const std::vector<double> &data
                                     )
{
    std::vector<double> PDF_x;
    std::vector<double> PDF_y;

    // CALCULATE CDF
    std::vector<double> CDF_x;
    std::vector<double> CDF_y;
    std::tie(
             CDF_x,
             CDF_y
             ) = distribution::CDF(
                                   data
                                   );

    // EVALUATE PDF BY FINITE DIFFERENCES

    // NOTE: The first point is estimated by a forwards finite-difference.
    //
    // NOTE: ... (This is preferred to a central finite-difference, by shifting the x point, in case the data is ordered, to keep the same x range, etc.)
    //
    double _PDF_x = CDF_x[0];
    double _PDF_y = (CDF_y[1] - CDF_y[0])/(CDF_x[1] - CDF_x[0]);

    PDF_x.push_back(_PDF_x);
    PDF_y.push_back(_PDF_y);

    // NOTE: ... The central points are estimated by central finite-differences.
    for( auto i = 1; i < (CDF_x.size() - 1); ++i )
    {
        _PDF_x = CDF_x[i];
        _PDF_y = (CDF_y[i+1] - CDF_y[i-1])/(CDF_x[i+1] - CDF_x[i-1]);

        PDF_x.push_back(_PDF_x);
        PDF_y.push_back(_PDF_y);
    }

    // NOTE: ... And the last point is estimated by a backwards finite-difference.
    _PDF_x = CDF_x.back();
    _PDF_y = (CDF_y[(CDF_y.size()-1)] - CDF_y[(CDF_y.size()-2)])/(CDF_x[(CDF_x.size()-1)] - CDF_x[(CDF_x.size()-2)]);

    PDF_x.push_back(_PDF_x);
    PDF_y.push_back(_PDF_y);

    // RETURN
    return std::make_pair(
                          PDF_x,
                          PDF_y
                          );
}
